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Abstract 



We consider the use of A'"-body simulations for studying the evolution of rich star 
clusters (i.e. globular clusters). The dynamical processes included in this study are re- 
stricted to gravitational (point-mass) interactions, the steady tidal field of a galaxy, and 
instantaneous mass loss resulting from stellar evolution. With evolution driven by these 
mechanisms, it is known that clusters fall roughly into two broad classes: those that dis- 
sipate promptly in the tidal field, as a result of mass loss, and those that survive long 
enough for their evolution to become dominated by two-body relaxation. 

The time scales of the processes we consider scale in different ways with the number 
of stars in the simulation, and the main aim of the paper is to suggest how the scaling of a 
simulation should be done so that the results are representative of the evolution of a "real" 
cluster. We investigate three different ways of scaling time. One of these is appropriate 
to the first type of cluster, i.e. those that dissipate rapidly, and similarly a second scaling 
is appropriate only to the second (relaxation-dominated) type. We also develop a hybrid 
scaling which is a satisfactory compromise for both types of cluster. Finally we present 
evidence that the widely used Fokker-Planck method produces models which are in good 
agreement with A'"-body models of those clusters which are relaxation dominated, at least 
for A'"-body models with several thousand particles, but that the Fokker-Planck models 
evolve too fast for clusters which dissipate promptly. 

Key words: gravitation - methods: numerical - celestial mechanics, stellar dynamics 
- globular clusters: general 
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1. Introduction 

This paper is concerned with techniques for studying the dynamical evolution of a glob- 
ular star cluster. This is a complicated problem: a complete description would be impossi- 
ble without a detailed understanding of the way in which stars interact non-gravitationally, 
and how the internal evolution of the stars interacts with their dynamical behaviour. This 
is illustrated by the emerging study of the behaviour of binaries in star clusters (Milone 
and Mermilliod 1996). Nevertheless there is no astrophysical problem which admits of a 
complete description, and much may be learned by the study of more tractable, idealised 
models, provided that the goal of understanding the behaviour of real clusters is always 
kept in mind. 

To a very good approximation the internal dynamics of a globular star cluster is an 
A^-body problem. The dominant forces which govern the internal motions (i.e. motions 
relative to the centre of mass of a cluster) are the gravitational forces of the stars themselves 
and the tidal field of the parent galaxy. In addition, the masses of the stars are constant for 
long periods of time. Nevertheless, stars lose significant amounts of mass from about the 
time at which they complete their main sequence evolution, and it is known from simplified 
models (e.g. Applegate 1986, Chernoff & Shapiro 1987) that the effects on the cluster can 
be devastating. It is fair to say that this is the principal known channel through which the 
internal evolution of the stars can influence the gross dynamics of a cluster. 

We have now stated the three main processes whose effects we explore in this paper: 

(i) mutual gravitational interactions of the stars; (ii) the tidal field of the galaxy; (iii) 
mass loss by stellar evolution. Clearly much has been omitted. For example, nothing 
has been said about the effects of primordial gas, whose early dissipation may well be 
sufficient to unbind a young cluster (e.g. Lada et al. 1984, Goodwin 1997). Nevertheless, 
as such studies show, the infiuence of this matter may well almost cease very early on in 
the evolution of a cluster, i.e. by the time the most massive stars have evolved. Another 
omission is the infiuence of primordial binaries (e.g. Hut et al. 1992). Their effect is more 
persistent, and indeed they may well be crucial for long periods in the later evolution of a 
cluster, probably determining such factors as the size of the core (Goodman & Hut 1989, 
Gao et al. 1991, Vesperini & Chernoff 1994). Nevertheless, up to the point of core collapse 
(the time scale for which is one of the features on which we concentrate in this paper), 
primordial binaries are thought to behave roughly as a species of rather massive stars, i.e. 
with a mass equal to the sum of the masses of the components (Heggie & Aarseth 1992). 
Up to this point, then, it is not expected that they have a very important effect on the 
gross dynamics. 

The simplified problem that we have stated can be tackled by more than one technique. 
At present the most popular is the Fokker-Planck model, which treats a cluster of stars 
rather like a gas, the gravitational encounters between the stars being modelled in much 
the same way as collisions in an ideal gas. This technique is idealised in a number of 
ways: often it is assumed (i) that the cluster and tidal field are spherically symmetric, 

(ii) that the distribution of stellar velocities is isotropic, (iii) that the evolutionary time 
scale is long compared with the time scale on which stars orbit within the cluster, and 
(iv) that small-angle scattering encounters are dominant. Techniques for dispensing with 
most of these assumptions are known (e.g. (i) Goodman 1983a, Einsel & Spurzem 1998; 
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(ii) Takahashi 1995; (iii) Spitzer & Hart 1971; (iv) Goodman 1983b), but they are not in 
routine use at present. 

Apart from Fokker-Planck techniques, another potential contender is the A^-body 
method (e.g. Aarseth 1985). This has no difficulty with any of issues (i)-(iv) raised above. 
Indeed there are important ways in which AT-body techniques become easier if the models 
are made more realistic. One example is the effect of introducing a spectrum of stellar 
masses. Multi-mass Fokker-Planck models are more time-consuming than models in which 
all stars have the same mass, whereas A^-body models are easier because the time scale for 
evolution becomes shorter. 

In general A^-body models are more free of simplifying assumptions than other tech- 
niques. Their great drawback is that it is still impractical to model systems with the 
correct A^, because of the computing time required. On a general-purpose supercomputer 
a 10,000-body model requires of order 2 months to pass core collapse (with equal masses; 
Spurzem & Aarseth 1996), while even the fastest special-purpose computer (the GRAPE-4 
in Tokyo) takes of order 3 months for about 32,000 particles (Makino 1996; the calculation 
was actually performed with only one quarter of the full configuration of the hardware). 
These are small simulations compared with a median globular cluster, which may con- 
tain of order 10^ stars. Therefore the main problem in studying the evolution of globular 
clusters with A^-body models is knowing how the results depend on A^. 

To some extent this question could be answered empirically, by computing simulations 
with different N and extrapolating. Theory, however, gives important information on 
the A"-dependence of known evolutionary processes such as two-body relaxation, and it 
is desirable that the extrapolation should be consistent with this. It follows that the 
application of A^-body techniques to the evolution of globular clusters is to some extent as 
dependent on theory as the Fokker-Planck model is. 

In this paper we have two aims. One aim is to show that the results of A^-body 
simulations can indeed be scaled reliably to systems of the size of globular clusters. We do 
so primarily by comparing simulations with different A^, the scaling being determined in 
accordance with theoretical considerations. Our second aim, however, is to demonstrate 
how the resulting evolution compares with that obtained with Fokker-Planck techniques 
incorporating the same physical mechanisms. 

To some extent both questions have already been studied by Fukushige and Heggie 
(1995), who compared the results of A"-body and Fokker-Planck calculations, and varied 
the number of particles in the simulations until the results were nearly independent of N. 
That paper imposed on itself two limitations, however. First, the paper was concerned 
only with clusters in which relaxation effects were unimportant: the models dissolved in 
the tidal field before there was any significant mass segregation. The second limitation, 
which derives from the first, was that the scaling of the models to real star clusters was 
uncomplicated, essentially because it was unnecessary to scale the results so that relaxation 
effects in the models would closely mimic those in a real cluster. 

In the present paper we remove these limitations, and extend the problem by consid- 
ering how to model those clusters in which relaxation is important, which in effect means 
those in which there is relatively little mass loss by stellar evolution, and/or the initial con- 
figuration is sufficiently centrally concentrated that relaxation effects can occur sufficiently 



4 



quickly. 

Here is the plan of the paper. In the next section we introduce the physical mechanisms 
which drive the evolution, summarise our choice of initial conditions, outline some of the 
hardware and software aspects of our simulations, and discuss the central issue of time 
scaling. Section three describes the results, which show that clusters whose evolution is 
driven by tidal overflow, two-body relaxation, and mass-loss from stellar evolution, may 
indeed be simulated with relatively modest values of A^. We also show that the results 
of simulations with a few thousand particles are remarkably consistent with those already 
found using Fokker-Planck techniques, at least for those clusters which are not so seriously 
aflFected by the third mechanism (i.e. mass loss by stellar evolution) that they disrupt in 
a time much less than a Hubble time. In the last section of the paper we summarise our 
conclusions, and also draw attention to some astrophysically interesting features of the 
evolution, such as the anisotropy of velocities. 

A companion paper (Vesperini & Heggie 1997) applies some of the techniques of the 
present paper to discuss the evolution of the mass function. Finally, the question of how to 
scale A'"-body models so that the statistics of stellar collisions in evolving globular clusters 
may be studied is the subject of a further paper (Aarseth & Heggie 1998). 

2. Technique 

2.1 Initial and Boundary Conditions 

As has been stated, one of our aims has been to compare the results of A^-body 
simulations with those of Fokker-Planck calculations. Fortunately there exists an excellent 
and instructive set of such calculations described by Chernoff & Weinberg (1990), and 
so we have adopted their initial and boundary conditions (with appropriate modifications 
described below), even where these differ from those that might now be considered more 
appropriate for the study of galactic globular clusters (e.g. the minimum mass). In one or 
two cases we have carried out a few A/"-body simulations to observe the effect of changing 
these parameters, but much more detailed results along these lines will be found in the 
paper by Vesperini & Heggie (1997). 

Chernoff & Weinberg studied models of clusters in circular orbit about a galaxy whose 
gravitational field is taken to be that of a point mass. The circular speed is Vg — 220km/s. 
In our A'"-body simulations we have adopted the corresponding equations 

X — 2(jjy — Zlu X = Fx 

y + 2Lox = Fy (1) 
z + uPz = Fz 

(e.g. Chandrasekhar 1942), where the coordinates of a star are taken relative to a rotating 
frame whose origin rotates about the galaxy in a circular orbit of angular velocity a;; the 
X- and y-axes point in the direction away from the galactic centre and in the direction 
of galactic rotation, respectively; and F is the acceleration due to other cluster members. 
Note that Fokker-Planck codes usually impose a cutoff, rather than a tidal field as in the 
A'"-body simulations, though the evidence is that this makes little difference to the rate at 
which mass is lost (Giersz & Heggie 1997). 
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For initial conditions ChernofF & Weinberg (1990) used King models (King 1966) with 
scaled central potential Wq = 1,3 and 7. These models have also been used in the present 
study, but with Wq =3,5 and 7. The reason for ignoring the models of lowest concentration 
is that Chernoff & Weinberg found that they always disrupted in much less than a Hubble 
time, a result that has already been confirmed (with substantial quantitative differences) 
with A'"-body models by Fukushige & Heggie (1995). Chernoff & Weinberg found that, 
depending on the initial mass function, clusters with Wq = 3 and 7 could survive until 
the present day, and we have added the case Wq = 5 to help reveal trends in the results. 
In planning our programme of work we regarded the results of Fukushige & Heggie as 
providing a definitive AT-body answer to the evolution of those clusters which dissolve 
quickly, without surviving long enough for the effects of relaxation to become important. 
Therefore we have concentrated on those models whose long-term fate is determined in an 
important way by relaxation effects, and was left open by Fukushige & Heggie. 

We have followed Chernoff & Weinberg in the choice of initial mass function, which 
is a power law dN oc m~°'dm in the range ^AMq < m < ISM©, with a = 1.5, 2.5 or 3.5. 
We have also followed their prescription for mass loss through stellar evolution, which is 
assumed to take place instantaneously at the end of main sequence evolution. The time 
and amount of mass lost are determined by linear interpolation in a table. At the time of 
its formation a remnant is assumed to have exactly the same velocity as its progenitor; in 
making this assumption we are again following Chernoff & Weinberg, though it is likely 
to be in serious error for neutron stars, only a small fraction of which would probably be 
retained (Drukier 1996). There are no primordial binary stars in our models, though these 
will be included in the next paper in this series (Aarseth & Heggie 1998). 

Let us consider next the number of parameters in this problem. There are two dimen- 
sionless parameters, Wq and a, which we retain, and two dimensional parameters, which 
are the initial mass of the cluster, M(0), and its galactocentric distance, Rg. As shown 
by Chernoff & Weinberg, clusters with the same relaxation time exhibit the same evolu- 
tion, if the assumptions of the orbit-averaged Fokker-Planck model are valid. Therefore 
the evolution should depend on M(0) and Rg only through the parameter F introduced 
by Chernoff & Weinberg, as this is a measure of the relaxation time. We follow Chernoff 
& Weinberg in presenting results for four distinct values of F, corresponding to the four 
"families" that they studied. In Sec. 3. 2 we check whether different models within the same 
family evolve similarly. 

2.2 Practical considerations 

In carrying out our A^-body models we have employed a version of the code NB0DY4 
(Aarseth 1985), which may be briefly described as follows. It is a direct summation code 
with a Hermite integrator (Makino & Aarseth 1992) and a binary hierarchy of time-steps 
(i.e. "block" time steps). It employs various mechanisms for the treatment of compact 
subsystems, including two-body and chain regularizations (Mikkola & Aarseth 1993). The 
inner binary of a hierarchical triple system is treated by a novel technique for rescaling of 
time (Mikkola & Aarseth 1996). The main feature of other advanced codes which NB0DY4 
lacks is the distinction between irregular and regular forces (Ahmad & Cohen 1973), i.e. 
those due to near neighbours and those due to the remainder of the system. 

The evaluation of forces and force derivatives (which are also required in the Hermite 
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integrator) has been carried out with the aid of a special-purpose computer. This is a 
single GRAPE-4 board containing 48 HARP processors (cf. Makino ct al. 1993), which 
also performs one or two other compute-intensive tasks, such as finding neighbours and 
potentials. The board connects to a standard Dec Alphastation 3000/700 via a HARP 
control board and interface. The performance of the board is summarised in Table 1; most 
of the results in this table are discussed in §3 below, but it also gives the CPU time for 
a variety of models. It can be seen that most models take less than one day, which is a 
reasonable target for a programme in which many models must be computed in order to 
explore an adequate range of parameters. 

It is interesting to compare the timing with that of the full-size GRAPE-4 in Tokyo 
(Makino et al. 1997). They timed equal-mass King models with Wq = 3 using a single 
cluster of nine processor boards, and we assume that our models with the same Wq and the 
steepest mass function {a = 3.5) are the most comparable. We find that our single board 
(with associated host) takes about 0.97 minutes for a single A^-body time unit (Heggie & 
Mathieu 1986) ii N = 4096, and this appears to be only about twice as long as on the nine 
boards in a single cluster in Tokyo. For N = 8192 the corresponding time for our single 
board is about 2.42 minutes, which exceeds the result of Makino et al. (on a cluster of nine 
boards) by a factor of a little more than two. In fact our timings are quite insensitive to 
the slope of the mass function. 

The fact that our system performs at almost half the speed of one whose peak speed 
is nine times higher can be interpreted in terms of efficiency, i.e the fraction of peak speed 
actually achieved. As shown by Makino et al. (their Fig. 17) the efficiency of their hardware 
is about 8% when N — 8192 (for a certain model). It may be expected that a smaller 
system, though slower, would be more efficient. Indeed, we estimate the efficiency of our 
hardware at about 25% for the same N, and this largely accounts for the relative timings. 

Now we turn to some practical aspects in the setting up of the initial A^-body model. 
Using the units introduced by King (1966), the desired King model was first constructed 
by numerical integration of Poisson's equation, which also yields the corresponding value 
of the tidal radius r^. Then an AT-body realisation of the model was constructed, also in 
King's units. This A"-body model was then scaled to standard A"-body units (Heggie Sz 
Mathieu 1986), in which the total mass of the cluster is 1, the constant of gravitation is 
1 and the initial energy is —1/4. Henceforth quantities expressed in these units will be 
distinguished by *. The scaling was carried out so that this A'-body model, regarded as 
an isolated model, is in virial equilibrium. The theoretical tidal radius of the model in 
King units was scaled in the same way, giving a value for the tidal radius of the A"-body 
model in A^-body units. Now in these units it follows from eqs.(l) that the theoretical 
tidal radius is also given by 

= (3cc;*2)-i/3^ (2) 

and so the value of uj* is determined. Note that the resulting initial model is not quite in 
virial equilibrium, because of the contribution to virial balance of the terms involving u 
in eqs.(l) (cf. Chandrasekhar 1942). This could be rectified (Heggie & Ramamani 1995), 
but the neglected contribution is no more than a few percent, and is not noticeable in 
comparison with statistical fluctuations in the virial ratio. 
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2.3 Scaling 



We now consider how the model so constructed is to be scaled in such a way that it 
represents the evolution of a "real" star cluster. As we shall see, the most delicate aspect 
of this question is how to scale the time from the A^-body units of the simulation to some 
astrophysical unit of time such as lO^yr. One of the principal roles of this scaling is that 
it determines the points in the simulation at which mass loss by stellar evolution takes 
place. In this subsection we first state two scalings which have an obvious physical basis, 
and then present a somewhat novel hybrid, which is designed to combine their advantages. 

2.3.1 Fixed Scalings 

The simplest choice is scaling by the crossing time, i.e. we write 



where tcr denotes the crossing time of the real cluster (in Myr, for example), and t*^ 
denotes that of the A^-body model (in A^-body units). This choice guarantees that the 
period of radial motion of a particle in the simulation scales to the period of radial motion 
of a star in the cluster. It may be shown readily that it also ensures that the angular 
velocity cu*, defined by eq.(2), scales correctly to the orbital angular velocity of the cluster, 
i.e. to a; = Vg/Rg. Thus eq.(3) is equivalent to assuming that t = t*Ut, where 



In general, scaling by the crossing time is appropriate at times when mass loss is important, 
because simple arguments show (Hills 1980) that the effect of impulsive mass loss, i.e. mass 
loss taking place on a time scale much smaller than the crossing time, is quantitatively 
very different from that of adiabatic mass loss. 

Scaling by the crossing time, eq.(3), does not ensure that the number of relaxation 
times that elapse in the two systems correspond correctly, because the relaxation time is 
an A^-dependent multiple of the crossing time, and the number of particles in a practical 
simulation is much smaller than that in the globular clusters of interest. A modified time 
scaling, which ensures that relaxation proceeds at the correct rate, is 



where trh is the half-mass relaxation time in the real cluster, and t*^ is the corresponding 
quantity for the A^-body model in A^-body units. (Any other measure of relaxation times, 
such as the central value, would be equally suitable.) 

2.3.2 Variable Scaling 

Scaling by the crossing time, i.e. eq.(3), was used by Fukushige & Heggie (1995) to 
simulate the early evolution of globular clusters, but this meant that these authors were 
unable to determine the long-term fate of the clusters which survived this early phase, i.e. 
those with steep initial mass functions and/or high initial concentration, for the evolution 




(3) 



(4) 




(5) 
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of such clusters becomes relaxation-dominated. It is important to be able to simulate such 
clusters, for if the initial conditions of the galactic globular clusters bear any resemblance 
to those of the models surveyed by Chernoff & Weinberg, it is these long-lived models 
which are relevant to the evolutionary history of the clusters which exist today. Ideally, 
therefore, one should devise a strategy which combines the advantages of both kinds of 
scaling, i.e. eq.(3) for a cluster which dissolves promptly, and eq.(5) for a long-lived cluster, 
at least during the long phase of its evolution which is dominated by relaxation. In order 
to decide which scaling to adopt, and when, we may compare the time scales for the two 
dominant processes, i.e. mass-loss and relaxation. 

The time scale for mass-loss by stellar evolution may be defined as 



tM = -M/{dM/dt), (6) 

where the derivative on the right indicates only mass loss arising from this mechanism 
(cf. eq.[8] and Fig.l). We denote the time scale of two-body relaxation by tr, and defer 
a precise definition until later in this section. Then our strategy suggests that we adopt 
eq.(5), i.e. scaling by the relaxation time, if Im > tr- We must be more careful, however, 
before concluding that we use eq.(3) if Im < tr, because this scaling does not correctly 
scale the relaxation time of the cluster to that of the model. Therefore a situation could 
arise in which t* < t\^, even though tr > tM- For this reason we introduce an intermediate 
regime in which we compromise as best we can, by adopting a scaling which ensures that 
t* = t\j. It provides a smooth transition from the scaling of eq.(3) to that of eq.(5) as the 
mass loss from stellar evolution gradually slows down. 

At this point we summarise the variable scaling that is adopted in the bulk of the 
simulations described in this paper. Because we must switch between one scaling and 
another we will not usually be able to express t as a simple multiple of t*, and instead 
must consider scalings in differential form. Thus we have 

t* 

if tM <tcr-^, 

''cr 

t* 

^f ^crT^ <tM < tr, (7) 
'^cr 

if tM ^ tr - 



( tr 



dt 
dt* 



t* 
ir 



< tt 



Now we must consider how this may be implemented in an A'"-body calculation. We 
assume that, at any time, there is a power-law mass function with N(m)dm oc m~°'dm 
for mi > m > mo. Then the total mass is M oc m^~" — Wg"" and if we assume (for this 
purpose only) that each star loses all its mass at the end of its evolution, then it is easy 
to show that 



nil 

^ ^ dlogt 
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This is easily computed if the current value of mi is known (i.e. the mass of a star which 
is just reaching the end of its evolution) and if we assume that a retains its initial value; 

dj log Tfl\ 

and the value — is easily computed from the table of lifetime versus mass. Fig.l 

alogt 

shows the result for several values of a. 

In eq.(7) the value of tu must be compared with tr-, and we must adopt a suitable 
definition of the relaxation time which allows us to determine whether the evolution of a 
cluster is dominated by relaxation or mass-loss. This is not straightforward. The evolution 
of the core of a high-concentration cluster may be dominated by relaxation even though 
mass-loss causes the half-mass radius to expand. The issue is further complicated by 
the mass-dependence of mass segregation, which depends in turn on the evolving mass 
function. Despite these complications we have adopted a simplistic criterion based on a 
result of Spitzer (1987, his eq.[4-l]) for values at the half-mass radius. We write it in the 
form 

= ^ (9) 
t,r llln(7iV)- ^ ' 

The current value of may be obtained from the current value of A^* and the initial values 
N{0) and iV*(0), where N{0) = M(0)/m(0), if we assume that N/N* is constant, and this 
is correct provided that the model correctly simulates the evolution of the real cluster. 

The remaining quantities in eq.(7) are easily calculated. Thus tcr = Utt^^, where 
Ut is given by cq.(4), and t*^ may be computed from the conventional definition t*^ = 
M*^/V(2|E* 1)3/2, in which E* is the energy of the model in A'"-body units. The ratio 
t*/tlj. is obtained in analogy with eq.(9). 

We close with a few practical details. First, we adopt a relatively low value of 7 ~ 0.01, 
based on the work of Giersz & Heggie (1997). Second, as N* can become very small, the 
factor ln7A'^* in these expressions was modified to ln(l -|- ^N*). When N* = 200, for 
example, the value of dt/dt* actually used is too large by a factor of about 1.6. Finally, 
for astrophysical reasons it was deemed sufficient to terminate each simulation by the time 
t = 20Gyr. 

3 Numerical Results 

3.1 Tests of the scaling strategy 

3.1.1 A'"-Dependence of Scaling 

In the previous section we have introduced three ways of scaling the time in a small 
A'^-body simulation to that of a typical globular cluster. We refer to these respectively 
as tc-scaling (equation 3), tr-scaling (equation 5) and variable scaling (equation 7). The 
most straightforward test of any of these is the comparison of A^-body simulations carried 
out with difi^erent values of N*. This is the purpose of the present subsection. In §3.2 we 
go on to carry out a comparison between the A'"-body results and those of Fokker-Planck 
models. 

Fig. 2 shows an example of the same model, specified by the information in the title, 
computed with different A^*, and using tj.-scaling (equation 5). The uppermost curve in 
each plot shows the tidal radius, and the remainder are 0.1%, 1%, 10% and 50% Lagrangian 
radii, as defined in the caption. From a comparison of the two plots it is gratifying to note 
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that core collapse ends at very nearly the same time (about lOGyr). Towards the end 
of the simulations, however, the larger of the two (Fig. 2b, A^* = 8192) loses mass at a 
relatively larger rate. This can be seen from the uppermost curve (the tidal radius), as 
rt oc M^/^. This is an example of a general trend which will be discussed in connection 
with Table 1 below. 

Table 1 shows a more extensive set of results, though these models were computed 
with variable scaling (equation 7), and each model is now summarised by only a small 
number of data. These data include the time at which the model either dissipated or was 
terminated, tmax, the time at the end of core collapse, tec, and the mass and half- mass 
radius (in AT-body units) at this time and at 15Gyr. Many of the models, however, do not 
reach core collapse, and of these models most also disappear before 15Gyr. 

For fixed Wq and a, the results from different A^* are in qualitative agreement. Quan- 
titatively, a systematic result emerges, which is that the larger model (A^* = 8192) tends 
to lose mass somewhat faster than the smaller model (cf. also Fig. 2). To quantify this we 
have determined the time (in Gyr) at which the value of M* in the smaller model takes the 
value at 15Gyr in the larger model. The resulting five values range from 15.3 to 17.4Gyr, 
with a mean of 16.6. This gives a quantitative estimate of about 10% for the order of 
magnitude of the discrepancy in the rate of loss of mass. 

3.1.2 Comparison of Different Scalings 

Though much was made in §2.3 about the desirability of variable scaling, in practice 
the differences between results produced by the different scaling methods are not usually 
very large, at least for initial particle numbers as large as N* = 8192. The basis for this 
assessment is the data in Table 2, which compares different types of scaling for the same 
clusters. Results for tc-scaling have been taken from Fukushige & Heggie (1995). 

Let us first concentrate on those clusters which dissipate promptly, either during or 
immediately following the period of heavy mass loss by stellar evolution. These are the 
three models in this table for which Fukushige & Heggie obtained a definitive result, viz. 
Wo = 3, a = 1.5 and 2.5, and Wq = 5, a = 1.5. For such clusters the dissolution time 
estimated by t^-scaling is systematically too large. The results from variable scaling are 
better, and we include here the case Wq — 3, a — 1.5, even though the lifetime given by 
variable scaling is considerably less than that of Fukushige & Heggie; their Fig. 6 suggests 
that statistical differences between models will lead to differences of order 0.03Gyr, and 
that the result stated in Table 2 (i.e. O.llGyr) is near the upper end of the range of values. 

Though variable scaling gives a better result for Wq = 3, q; = 2.5 than t^-scaling, it 
is still poor. Inspection of Fig. 7 in Fukushige & Heggie (1995) suggests that their result 
is quite robust in this case. Though data is not given in Table 2, we have found that, for 
smaller particle numbers, the disagreement is still larger, as expected. 

The model that we have been discussing {Wq = 3, q; = 2.5) is expected to be diffi- 
cult to simulate, because it lies close to the border between those clusters which dissipate 
promptly (without significant relaxation) and those that can survive for several Gyr. An- 
other problematic case of this kind is Wq = 7, a = 1.5. Fukushige & Heggie stopped their 
calculation at 1.6Gyr, when the mass of this cluster was approximately M* — 0.05 (their 
Fig.l), whereas our model with variable scaling dissolved at 1.2Gyr. It seems likely that 
the model of Fukushige & Heggie, despite the softening they used, is significantly affected 
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by two-body relaxation. This is the longest of their simulations (to t* = 3000), and their 
Fig. 8 shows that even a shorter run (to t* = 1000, though with different model parameters) 
may be significantly affected by two-body relaxation. If so, the resulting core evolution 
could make the cluster too robust to the effects of mass loss, and so prolong its life. 

For the remaining five sets of models presented in Table 2 (which, as it happens, are 
those that survive to 15Gyr) the results are more straightforward. For these models the 
work of Fukushige & Heggic gives only a lower limit, and can be ignored. The results from 
variable- and tj,-scaling are quite consistent, the only noticeable general trend being that 
the models with t^-scaling evolve more slowly. The one exception to this trend is the core 
collapse time for Wq = 7, a = 2.5, but this quantity is known to be subject to considerable 
statistical fluctuation (Spurzem & Aarseth 1996). 

3.2 Comparison with Fokker-Planck Results 

In this subsection we consider the second main purpose of the investigation, which is 
to compare A/^-body results with Fokker-Planck results for the same cluster parameters, 
partly as a consistency check, and partly as a check on the validity of the Fokker-Planck 
model, which is likely to remain a standard tool for the investigation of the dynamical 
evolution of globular clusters for some time to come. 

The basic information we present is given in Tables 3 and 4, which compare a number 
of data for our A'^-body models and a subset of the Fokker-Planck models of Chernoff 
& Weinberg (1990): we omit their models in which Wq = 1 initially. Note also that our 
models with Wq = 5 are omitted here because this was not one of the values which Chernoff 
& Weinberg included. All AT-body models use variable scaling (eq.(7)). 

A factor to be borne in mind in interpreting the Fokker-Planck data is that the time 
tmax does not correspond to the time at which the model loses all its mass. As Chernoff & 
Weinberg explain, there comes a point at which the method fails to find a new structure 
corresponding to the evolving distribution function, and this is interpreted as indicating 
the imminent rapid dissolution of the model. 

The caption in Table 3 (and those in Table 4, to be discussed below) refers to one of 
the "families" adopted by Chernoff & Weinberg (cf.§2.1). For the orbit-averaged technique 
used by Chernoff & Weinberg, the evolution is the same for all members of one family, 
given the initial value of Wq and the initial mass function. Thus a low-mass cluster at large 
galactocentric radius will have the same predicted evolution as a cluster of high mass at 
an appropriate small radius. Comparison of the models in Table 3 illustrates this scaling 
in the A^-body data, along with a comparison with the Fokker-Planck data. In each triplet 
of models the two A^-body models belong to the same family of Chernoff & Weinberg 
(family 1), but at different galactocentric radii. To avoid complications with the residual 
A'"-dependence of the results of the A'"-body models (cf.§3.1.1), all models in this Table 
share the same initial value of N*, taken to be 4096. The data of Table 3 shows yet again 
that the models divide into two types: those that dissolve in much less than a Hubble 
time, and those in which the A^-body models survive to at least 20Gyr. 

Let us consider first the rapidly dissolving models, i.e. those with {Wq, a) = (3, 1.5), 
(3, 2.5) and (7, 1.5). Though the A^-body models indicate systematically longer lifetimes at 
larger radii, the evidence is in fact not significant. We have already argued (in discussing 
Table 2) that the difference in the lifetimes of the two models with Wq = 3, q; = 1.5, 
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would probably be within statistical fluctuations even for simulations with N* = 8192. A 
similar argument probably applies to the A^-body models with Wq = 7, a = 1.5, though 
for somewhat more subtle reasons. These models rapidly lose most of their mass, retaining 
only about 20% by 0.5Gyr. The subsequent evolution is therefore subject to the sizeable 
fluctuations between different A-body models with rather small N*. 

Of the three models which dissolve rapidly, the Fokker-Planck and A-body results are 
consistent only for the model we have just discussed, i.e. Wq = 7, a = 1.5. The other 
two models {Wq = 3, a = 1.5 and 2.5) confirm the finding of Fukushige & Heggie (1995) 
that the Fokker-Planck model underestimates the true lifetime by a factor of order 10 in 
some cases. As they suggested, a plausible reason for this is that the Fokker-Planck models 
necessarily assume that Im, the time scale on which the cluster loses mass, much exceeds 
tcr, and that the lifetime is underestimated when this assumption is not valid. 

Now we turn to the more straightforward comparison of the three models in Table 3 
which survive for at least a Hubble time. The A-body models illustrate rather well the 
fact that the evolution is similar for clusters within the same family but at different Rg. 
The results also agree rather well with the Fokker-Planck data. Though it might seem that 
the Fokker-Planck models evolve somewhat faster than these A-body models, it must be 
recalled (cf. discussion of Table 1) that larger A-body models would also do so. 

Table 4 is the main vehicle for examining the relation between the A-body data and 
the Fokker-Planck data. Here all models are placed at about 4kpc, as the evidence of 
Table 3 suggests that the results will be independent of Rg within a given family, at least 
for long-lasting models. In fact we have omitted those models that dissolve rapidly, as an 
adequate comparison has been provided by Fukushige & Heggie. 

The data on core collapse in this Table is rather meagre and adds little to what was 
learned from Table 3, except that the models for Family 1 have larger A*. It confirms 
that there appear to be no systematic differences in either the time or the total mass at 
core collapse. 

Now let us consider the mass at 15Gyr. The absence of this data in two of the 
Fokker-Planck models of Family 1 is due to the fact that Chernoff & Weinberg stopped 
their calculations at the time of core collapse, which occurs before 15Gyr for these models. 
For the one remaining model in this family, and all models in the other families, there 
is evidence of some systematic differences in the mass at 15Gyr. In particular, the dis- 
agreement is most pronounced for the models with Wq — 3 and a = 3.5, especially for the 
first two families, i.e. those in which two-body relaxation is fastest. The Fokker-Planck 
models yield substantially smaller masses at 15Gyr. If we assume that the mass decreases 
approximately linearly with time, the differences correspond to differences in the lifetime 
which are in the range of 13-23%. The comparisons in Table 2, however, suggest that 
the disagreement is well within the residual A*-depcndcnce of the A-body results, and 
has the expected sign: the mass at 15Gyr decreases slightly with increasing A*, and the 
Fokker-Planck results are also smaller than the A-body results. These remarks apply to 
the models with Wo = 3 and a = 3.5; for other models the discrepancy in the mass at 
15Gyr is much more modest. 
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4. Discussion and Conclusions 
4.1 Anisotropy 

Though the main purposes of this paper are to investigate the scahng of AT-body 
models and to compare results with Fokker-Planck data, it is a pity to leave the N- 
body models without drawing attention to some other interesting results. One of these is 
the evolution of the mass function and the proportion of degenerate stars, and we refer 
interested readers to Vesperini & Heggie (1997), where detailed and comprehensive results 
are presented (though with a more realistic lower limit to the initial mass function, and 
disk shocking). 

Here we concentrate on the anisotropy of the models, which Chernoff & Weinberg 
did not address, as they used an isotropic Fokker-Planck model. For the purposes of this 
discussion we adopt the usual parameter (3—1 — where is half of the square 

of the transverse component of velocity of one star and Vr is the radial component, i.e. in 
the direction away from the density centre. Only stars inside the tidal radius are included 
in the averaging, which is not mass-weighted. Initially there is a short period in which 
/3 > 0, and this is attributable to the radial outward drift of stars during the early phase 
of rapid mass loss. Thereafter the sense of the anisotropy reverses, i.e. (3 < 0, and there is 
a slight predominance of transverse velocities. Nevertheless the anisotropy is small, values 
of order /? ~ —0.1 being typical. 

In an isolated model, (3 would become positive (e.g. Giersz & Heggie 1996), and the 
negative values in our model are attributable to the tidal field. This has two effects: in the 
first place it alters the angular momentum of the stars, as it is not spherically symmetric 
(cf.eq.[l]); and, secondly, it preferentially removes stars on nearly radial orbits, as these 
are most likely to escape. We attribute the negative values of (3 to the second effect 
mainly, as similar negative values of /3 are observed in models (not otherwise discussed in 
this paper) in which the tidal field of eq.(l) is replaced by a spherical cutoff, as in most 
implementations of the Fokker-Planck model. 

The smallness of the values oi \/3\ helps to justify the use of an isotropic model such 
as the Fokker-Planck model of Chernoff & Weinberg. On the other hand, our result also 
has a variety of implications for the interpretation of observational data. There are a few 
clusters for which careful study of proper motions requires the presence of significant radial 
anisotropy, i.e. (3 > (e.g. Cudworth 1979 for the cluster M3). There is a larger sample 
of clusters where careful dynamical modelling (using multi-mass anisotropic King-Michie 
models) again requires the presence of radial anisotropy (e.g. Meylan & Mayor 1991). Here 
the requirement of substantial anisotropy is less direct, as the models which are fitted are 
constrained by theoretical considerations (e.g. the choice of distribution function), and 
in principle it is possible that an incorrect choice of distribution function is compensated 
by the anisotropy when fitting to the observational data. Nevertheless, if the evidence for 
anisotropy is accepted, the contradiction with the anisotropy found in our models may 
only indicate that the initial conditions adopted in our survey differ too widely from those 
appropriate to the observed clusters. For example, all our models are initially tidally 
limited: the edge of the model is placed initially at the tidal radius. It is possible that 
initially much more compact models would develop strongly positive anisotropy. They 
would behave much more like isolated models until, after sufficient expansion in the post- 
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collapse phase of evolution, their development might more closely resemble that of the 
models we have studied. Interestingly, it was concluded by Phinney (1993), on quite 
different grounds, that only an initially very condensed cluster would evolve into an object 
resembling the cluster M15 at the present day. 

4.2 Conclusions 

Now we close with a summary of the main conclusions of this paper. It has been 
concerned with the dynamical evolution, over a time up to 20Gyr, of a set of models for 
the globular star clusters of the Galaxy. The evolutionary mechanisms considered here 
include all gravitational interactions between the members of a cluster, the tidal field of 
the Galaxy, and instantaneous mass loss at the end of main sequence evolution. We do 
not include the effects of finite stellar radii and collisions, or primordial binaries. The 
initial mass function is a power law between masses of 15 and OAMq, with index a. The 
initial model is a King model, whose concentration is specified by the dimensionless central 
potential Wq. 

In qualitative agreement with the results of Chernoff & Weinberg (1990) we find that 
the fate of the clusters divides them into two classes: those that disrupt relatively rapidly 
(in at most one or two Gyr) as a result of mass loss by stellar evolution and tidal stripping, 
without being significantly affected by two-body relaxation; and those that survive this 
initial phase of heavy mass loss and evolve towards core collapse. For clusters of initial 
mass about 1.5 x lO^M© and galactocentric radius of about 4kpc, clusters which disrupt 
are those with Wq = 3 and a ^ 2.5, and those with Wq = 5 or 7 and a ^ 1.5. In general 
those with larger Wq and/or larger a enter the phase in which their evolution is dominated 
by relaxation. 

For those clusters which disrupt promptly, a comparison between A'^-body and Fokker- 
Planck models was carried out by Fukushige & Hcggic (1995). In this paper we have 
concentrated on the longer-lived systems which become relaxation-dominated. For some 
parameter values the Fokker-Planck and A^-body results agree remarkably well; the time 
of core collapse of models with Wq = 7 and a ^ 2.5 agree as well as could be expected, 
bearing in mind the fact that there will be different collapse times even for different A^-body 
realisations with the same initial parameters. More systematic differences, of order 10%, 
arise for models which lie closer to the borderline between early disruption and relaxation- 
dominated evolution, e.g. the model with Wq — 3 and a — 3.5. For these models several 
lines of evidence (e.g. the mass at the present day, taken here as 15Gyr) imply that the 
Fokker-Planck models evolve somewhat more rapidly than the A^-body models. The sense 
of the disagreement is consistent with the conclusion of Fukushige & Heggie for the clusters 
which disrupt quickly, though the size of the disagreement is much smaller for the clusters 
which do not dissolve rapidly. Indeed the residual A^*-dependence of the rate of escape 
from the A'"-body models may explain part or all of the discrepancy. 

Besides a comparison with Fokker-Planck models, the other main purpose of this 
paper was an investigation of the way in which the evolution of small A^-body models may 
be scaled (in time) to that of globular clusters, where the number of stars is much larger. 
We have considered three possible scalings: (i) a scaling which ensures that the crossing 
time of the A'"-body model scales correctly to that of a real cluster; (ii) a scaling which 
ensures that the model relaxes at the same rate as the real cluster; and (iii) a hybrid, or 
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variable, scaling, which progresses from the first type of scaling to the second as the mass 
loss from stellar evolution slows down. 

The first type of scaling was applied by Fukushige & Heggie. It successfully allows the 
scaling of A^-body data for clusters which dissipate promptly without significant relaxation. 
The model relaxes too quickly, however, and so the scaling of the simulation is no longer 
valid as soon as relaxation in the model becomes significant. The second type of scaling 
produces results which appear to be correct for clusters which survive the initial phase of 
mass loss and enter the relaxation dominated phase of evolution. Unless the number of 
particles in the simulation is quite large, however, it can produce misleading results for 
those clusters which dissipate promptly. 

The third (variable) scaling was devised to incorporate the advantages of both fixed 
types of scaling. It produces qualitatively correct results for both kinds of clusters, i.e. 
those that dissipate rapidly and those that become relaxation-dominated. For clusters of 
the latter type its quantitative results are indistinguishable from those obtained by scaling 
of type (ii) . It is the most successful scaling method for clusters on the borderline between 
rapid destruction and long-term survival, e.g. King models with Wq = 7 and a = 1.5 (for 
our chosen initial range of stellar masses) . 

We recommend variable scaling as a means for using small A^-body simulations to 
model the evolution of globular star clusters. Its advantage is that it can be used (with 
some care) to model the evolution of clusters which dissipate promptly as well as those 
with long lifetimes. Though this might not seem an important property for the galactic 
globular cluster system, it is relevant if A^-body models are to be used for studying the 
evolution of other globular cluster systems, as in the Magellanic Clouds, which contain 
both young and old clusters. 

The results obtained in this paper are based on models which incorporate only a 
subset of the many dynamical and evolutionary processes which are needed in a full N- 
body description of a globular cluster (e.g. Aarscth 1996), but the problem of how to 
scale models cannot be avoided as long as the number of particles in a feasible A^-body 
simulation is much smaller than in a real globular star cluster. In a further paper (Aarseth 
& Heggie 1998) we shall consider how the radii of the stars may be scaled in order to 
simulate stellar collisions. 
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Table 1 



A'"-Dependence of Models with Variable Scaling 
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0.617 
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0.435 


0.94 



Notes 

1. An asterisk in the top line denotes a quantity in A'"-body units in an A'"-body model. 
The values of Wq, a and N* are initial values. Quantities at core collapse arc denoted by 
the subscript cc, and at 15Gyr by a subscript 15. The half-mass radius (of all particles 
inside the tidal radius) is denoted by r/j. 

2. The simulations end at tmax, where this value is less than 20, the calculation stopped 
with less than 10 bound particles, except that a plus indicates that it ended with a software 
or hardware error. 

3. In all these models the initial mass is M = 1.49 x 10^ Mq. 

4. This model was mainly computed without the HARP board. 
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Table 2 



Comparison of Different Time Scalings 
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Notes 

1. Most of the notation is defined in the notes to Table 1. 

2. Under the heading "Scaling", tc denotes scaling by the crossing time, eq.(3); tr denotes 
scaling by the relaxation time, eq.(5); and "var." denotes variable scaling, eq.(7). Results 
for ic-scaling are taken from Fukushige & Heggie (1995), and for several models the lifetime 
is only a lower limit. 

3. In all these models the initial mass is M = 1.49 x lO^M©. AU runs had N* = 8192 
particles initially. 
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Table 3 



Comparison of Fokker-Planck and A'"-body results for Familyl 
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7 


2.5 


FP 


9.6 




0.26 




7 


3.5 


4.1 


10.4 


20 


0.606 


0.476 


7 


3.5 


10.4 


10.5 


20 


0.62 


0.494 


7 


3.5 


FP 


10.5 




0.57 





Notes 

1. Most of the notation is defined in the notes to Table 1. 

2. Under the heading "Model", a number indicates results of an A^-body model (with 
iV* = 4096 initially) at the stated galacto centric distance (in kpc). "FP" denotes a Fokker- 
Planck model of Chernoff & Weinberg (1990). Their models were stopped at the point of 
dissolution of the cluster (in which case we give the corresponding time as tmax) ot at the 
end of core collapse (at time tec), whichever is earlier. 

3. The initial cluster mass is 5.45 x lO^M© at around lOkpc, and 1.49 x IO^Mq at around 
4kpc. All A'"-body models use variable scaling (eq.[7]). 
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Table 4 



Comparison with Fokker-Planck Results 



Wo 


a 


Rg (kpc) 


Model 


tec (Gyr) 










Family 1: 


M(0) = 


1.49 X lO^M© 






3 


3.5 


4.1 


8192 


> 18.5 




0.355 


3 


3.5 


4.1 


FP 


21.5 


0.078 


0.23 


7 


1.5 


3.7 


8192 


(1.2D) 





— 


7 


1.5 


3.7 


FP 


(l.OD) 


(0.022) 


— 


7 


2.5 


4.0 


8192 


11 


0.210 


0.089 


7 


2.5 


4.0 


FP 


9.6 


0.26 


- 


7 


3.5 


4.1 


8192 


9.2 


0.617 


0.435 


7 


3.5 


4.1 


FP 


10.5 


0.57 


- 






Family 2: 


M(0) = 


4.33 X lO^M© 






3 


3.5 


4.1 


4096 


> 20 




0.589 


3 


3.5 


4.1 


FP 


44.4 


0.035 


0.48 


7 


1.5 


3.7 


4096 


(2.0D) 





— 


7 


1.5 


3.7 


FP 


(3.0D) 


(0.0033) 


- 


7 


2.5 


4.0 


4096 


> 20 


— 


0.331 


7 


2.5 


4.0 


FP 


22.5 


0.26 


0.35 


7 


3.5 


4.1 


4096 


> 20 





0.698 


7 


3.5 


4.1 


FP 


31.1 


0.51 


0.68 






Family 3: 


M(0) = 


7.65 X lO^M© 






3 


3.5 


4.1 


4096 


> 20 




0.633 


3 


3.5 


4.1 


FP 


(42.3D) 


(0.085) 


0.58 


7 


1.5 


3.7 


4096 


(3.1D) 





— 


7 


1.5 


3.7 


FP 


(4.2D) 


(0.0080) 


- 


7 


2.5 


4.0 


4096 


> 20 


— 


0.374 


7 


2.5 


4.0 


FP 


35.5 


0.26 


0.40 


7 


3.5 


4.1 


4096 


> 20 


— 


0.741 


7 


3.5 


4.1 


FP 


51.3 


0.48 


0.71 






Family 4: 


M(0) = 


2.17 X lO^M© 






q 
o 


o.o 


A 1 


^uyo 






n fiQ9 


3 


3.5 


4.1 


FP 


(43.5D) 


(0.28) 


0.64 


7 


1.5 


3.7 


4096 


(5.0D) 






7 


1.5 


3.7 


FP 


(5.9D) 


(0.023) 




7 


2.5 


4.0 


4096 


> 20 




0.443 


7 


2.5 


4.0 


FP 


83.1 


0.25 


0.45 


7 


3.5 


4.1 


4096 


> 20 




0.783 


7 


3.5 


4.1 


FP 


131.3 


0.49 


0.77 
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Notes 

1. Most of the notation is explained in the notes to the previous tables. 

2. A number in the column "Model" indicates that the results are for an A^-body model 
with the initial value of N* as stated. 

3. Though the time of core collapse, tec, is given as "> 20" for those models which stopped 
at 20Gyr without collapsing, it is not possible to state whether these models would reach 
core collapse before dissipating. Where the value in this column is enclosed in brackets, 
with a "D", the model dissipated at the time stated, without previously undergoing core 
collapse; for Fokker-Planck models the mass at this time is entered in parentheses in the 
next column. 
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Figure Captions 



Fig.l Time scale of mass loss (eq.(8)) for mass functions with minimum mass OAMq and 
mass function index a — 1.5, 2.5 and 3.5 (where Salpeter is 2.35). The piecewise linear 
nature of the curves arises from low-order interpolation in the table of evolution time as a 
function of initial stellar mass (Chernoff & Weinberg 1990). 

Fig. 2 Comparison of results for the same initial and boundary conditions, but with different 
initial values of N*: (a) N* = 4096, (b) A^* = 8192. The tidal radius and four Lagrangian 
radii (corresponding to the innermost 0.1%, 1%, 10% and 50% of the mass within the tidal 
radius, measured from the density centre) are plotted against time. 
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